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ABSTRACT 

Roche tomography is a technique used for imaging the Roche-lobe filling secondary 
stars in cataclysmic variables (CVs). In order to interpret Roche tomograms correctly, 
one must determine whether features in the reconstruction are real, or due to statistical 
or systematic errors. We explore the effects of systematic errors using reconstructions of 
simulated datasets and show that systematic errors result in characteristic distortions 
of the final reconstructions that can be identified and corrected. In addition, we present 
a new method of estimating statistical errors on tomographic reconstructions using a 
Monte-Carlo bootstrapping algorithm and show this method to be much more reliable 
than Monte-Carlo methods which 'jiggle' the data points in accordance with the size 
of their error bars. 
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1 INTRODUCTION 

The secondary, Roche-lobe filling stars in CVs are key to 
our understanding of the origin, evolution and behaviour 
of this class of interacting binary. To best study the sec- 
ondary stars in CVs we would ideally like direct images of 
the stellar surface. This is currently impossible, however, as 
typical CV secondary stars have radii of 400 000 km and 
distances of 200 parsecs, which means that to detect a fea- 
ture covering 20 per cent of the star's surface requires a 
resolution of approximately 1 micro-arcsecond, 10 000 times 
greater than the d iffraction limited resolutio n of the world's 
largest telescopes. |R,utten fc Dhillon (1994) described a way 
around this problem using an indirect imaging technique 
called Roche tomography which uses phase-resolved spectra 
to reconstruct the line intensity distribution on the surface 
of the secondary star. 

Obtaining surface images of the secondary star in CVs 
has far-reaching implications. For example, a knowledge of 
the irradiation pattern on the inner hemisphere of the sec- 
ondary star in CVs is essential if one is to calculate stel- 
lar masses a ccurately enough to tes t binary star evolution 



models (see Smith fc Dhillon 1998). Furthermore, the ir- 



radiation pattern provides information on the geom etry of 



the a ccreting structures around the white dwarf (see Smith 
1995). Perhaps even more importantly, surface images of 



novae, cataclysmic variables - stars: late-type - stars: 



CV secondaries can be used to study the solar-stellar con- 
nection. It is well known that magnetic activity in isolated 
lower-main seque nce stars inc reases with decreasing rota- 
tion period (e.g. Rutten 1987). The most rapidly rotating 
isolated stars of this type have rotation periods of ~ 8 
hours, much slower than the synchronously rotating sec- 
ondary stars found in most CVs. One would therefore ex- 
pect CVs to show even higher levels of magnetic activity. 
There is a great deal of indirect evidence for magnetic ac- 
tivity in CVs - magnetic activity cycles have been invoked 
to explain variations in the orbital periods, mean bright- 
nesses and m ean intervals between outbursts in CVs (see 
Warner 1995| ). The magnetic field of the secondary star is 
also believed to play a crucial role in angular momentum 
loss via magnetic braking in longer-period CVs, enabling 
CVs to transfer mass and evolve to shorter periods. One of 
the observable consequences of magnetic activity are star- 
spots, and their number, size, distribution and variability, 
as deduced from surface images of CV secondaries, would 
provide critical tests of stellar dynamo models in a hitherto 
untested period regime. 

However, determining whether features in Roche to- 
mography reconstructions are real or not is problematic as 
the resulting maps are prone to both systematic errors, due 
to errors in the assumptions underlying the technique, and 
statistical errors, due to measurement errors on the observed 
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data points. It is therefore essential to quantify the effects of line profile models ( Unruh fc Collier Cameron 1995|) , chro- 
these e rrors on Roche tomograms in order to properly assess 
the rea 



ity of any features present. In this paper we study 



mospheric emission ( Dhruh fe Collier Cameron 1997j Bruls 
Solanki & Schussler 1998]), gravity darkening and differential 



the artefacts produced by errors in the input parameters 
using reconstructions of simulated datasets, and present a 
techniq ue for estimating the statistical errors on the recon- 



rotation ( patzes*^t**aTl996 ). In this section we explore the 
effects of systematic errors on Roche tomograms of C Vs us- 
ing simulated reconstructions in a similar manner to (Marsh 



structic iis. Subsequent papers will present the application fc Home (1988)| in their treatment of systematic errors 



of Roche tomography to real datasets. 



Doppler tomograms of CV accretion discs. 



3.1 The test image 



2 ROCHE TOMOGRAPHY 



Roche tomog raphy (Rutten & Dhillon 1994; Rutten & 



Dhillon 1996) is analogous to the Doppler imagi ng tech- 
nique used to map rapidly rotating sin gle stars (e.g. Vogt & 
Pernod 1983|), detached binary stars ( Ramscycr, Hatzes & 
Jablonski 1995) and contact binaries (Maceroni et al. 1994[). 



However, in contrast to the Doppler imaging technique, the 
continuum is ignored in Roche tomography because of the 
unknown and variable contribution of the accretion regions 
to the spectrum, which forces Roche tomography to map ab- 
solute line fluxes. Apart from this, the principles only differ 
in detail. 

In Roche tomography we assume values for the binary 
parameters and that the secondary is Roche-lobe filling, 
locked in synchronous rotation and has a circularised orbit. 
We then model the secondary as a series of quadrilateral 
tiles of approximately equal area lying on the critical Roche 
surface. Each tile or surface element is assigned a copy of the 
local (intrinsic) specific intensity profile convolved with the 
instrumental resolution. These profiles are scaled to take into 
account the projected area, limb darkening and obscuration 
and then Doppler shifted according to the radial velocity 
of the surface element at a particular phase. Summing up 
the contributions from each element gives the rotationally 
broadened profile at that particular orbital phase. 

In order to reconstruct a surface image of the secondary, 
the 'inverse' of this process must be carried out. We do this 
by iteratively varying the strengths of the profile contributed 
by each element assuming that the shape of the intrinsic line 
profile does not change. This iterative procedure is carried 
out (assuming no contribution from the accreting regions) 
until a satisfactory fit to the observed data, defined by the 
X 2 statistic (i.e. \ 2 ~ 1)) is achieved. As there are a large 
number of different intensity distributions that all predict 
trailed spectra consistent with the observed one, an addi- 
tional constraint is required. This is found by selecting the 
map of maximum entropy with respect to a suitable default 
map and is perfor med using the maximum entropy algo- 
rithm developed by Skilling fc Bryan (1984) . Further details 



on the Roche tom ography process can be found in 
fc Watson ( 2000)|. 



Dhillon 



3 SYSTEMATIC ERRORS AND ARTEFACTS 

There have been a number of studies exploring the effects of 
systematic errors on Doppler images of single stars. For ex- 
ample, the significance of so-called polar spots (large, long- 
lived, high latitude spots straddling the polar regions) has 
been tested for robustness against errors such as incorrect 



All our simulations are based on a single test image (Image 
A, figure |l|). This image mimics the effect of irradiation of 
the secondary star by the compact primary, taking into ac- 
count the distance to the companion and the incidence angle 
at the surface of the star. The intensity of the side of the 
star that is not irradiated is set to the intensity at the termi- 
nator. Also, shielding by (for example) an accretion stream 
is mimicked by a narrow band stretching along the equator 
from the inner Lagrangian (Li) point to the terminator and 
covering 3.6% of the total surface area. The intensity of this 
band is equal to the intensity of the non-irradiated side of 
the secondary. 

In addition, an array of ten completely dark spots are 
distributed around the star. Four are located on the trail- 
ing hemisphere centered at latitudes of approximately +50° , 
+25°, 0° and -25°, all with different longitudes. Another 
four are located on the leading hemisphere at latitudes of 
+50°, +25°, 0° and -36°, all with the same longitude. The 
two final spots are located at the north pole (+90°) and the 
Li point. Each spot covers 0.7% of the total surface area 
of the secondary, with the exception of the Li spot which 
covers 1%. 

Finally, the geometry of the Roche-lobe itself is given 
by the binary parameters Mi = 1.0 Mq, M2 = 0.5 Mq and 
an orbital period of 2 hours. Unless otherwise stated the 
orbital inclination used was 60°. In total, the image con- 
sists of 3008 elements which, with these binary parameters, 
gives a maximum radial velocity difference between any two 
neighbouring elements of 18.1 kms -1 . 

3.2 Synthetic datasets and fits 

Unless otherwise stated, a perfect synthetic trailed spectrum 
was computed assuming an instrumental resolution of 20 
kms -1 and a negligible intrinsic line width. The whole or- 
bit was sampled over 50 evenly spaced phased bins, with a 
velocity bin- width of 10 kms -1 and a negligible exposure 
time. No contribution from an accretion disc or other part 
of the system was considered. A sample trailed spectrum is 
shown in figure ^ and the best fit to this dataset is displayed 
in Image B, figure |lj 

In each case the initial guess at the fit was one of uni- 
form intensity distribu tion and a uniform default map was 
used except in section 3.14. Due to the nature of the tests, 
the fits could not proceed to the same final \ 2 on each oc- 
casion. Thus the optimal fit was determined by minimising 
X 2 until the spot features were resolved. 

The reconstructions are plotted using a linear grey- 
scale. A value of 100 corresponds to the maximum intensity 
on the original test image (Image A, figure [l]). Any inten- 
sity value ~ 22% greater than the maximum intensity of 
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Figure 1. The effects of systematic errors on the test image. 
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Image J. The 'mirroring' effect 
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Figure 1 (continued). The effects of systematic errors on the test image. 
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velocity (km/s) 

Figure 2. Synthesized dataset produced by the test image. 

the test image is flagged as white. This is because some 
artefacts generate 'hot' elements with high intensities which 
would otherwise dominate the grey-scale plots. 

We are now ready to explore the effects of systematic 
errors on the Roche tomography reconstructions. In each 
case a fake dataset was generated that was either corrupted 
in some form or was fitted assuming an incorrect parameter. 



3.3 Systemic velocity 7 

The systemic velocity is the radial velocity of the centre-of- 
mass of the binary system and represents the mean shift of 
the line from its rest wavelength. In cases where the wrong 
value for the systemic velocity is assumed, the reconstruction 
process will have difficulty in fitting the line profiles. As 
a result, artefacts will be generated in order to allow the 
model to fit the profiles as well as possible. Due to the high 
quality of the synthetic datasets used, one might expect a 
small error in the systemic velocity to produce noticeable 
artefacts. Image C (figure shows the effects of ignoring a 
systemic velocity of +5 kms -1 present in the data. Although 
the spot features and, to some extent, the irradiated inner 
face of the secondary are still visible in the reconstruction, 
dark bands framed by bright stripes are formed around the 
equatorial regions. 

These equatorial stripes can be explained as follows. 
When fitting the trailed spectrum the data appears to be 
'skewed' to one side of the radial velocity axis and therefore 
the model cannot possibly fit the edges of the line profiles. 
At the lower (or most negative) radial velocities expected by 
the model there is no data, the data that should have been 
present at these velocities having been shifted to a higher 
radial velocity. Thus the model attempts to fit zero intensity 
to the regions of the star that contribute to the line profile 
at these points, resulting in a dark equatorial band. 

At the other extreme, the higher (or most positive) ra- 
dial velocities in the line profiles are shifted beyond the ve- 
locities expected by the model. Thus the intensities found 
at the highest radial velocities are far greater than expected 
by the model and the opposite scenario emerges. The model 
now attempts to fit high intensities to the regions that con- 
tribute to the line profiles at these points. As only the equa- 
torial regions ever contribute to the edges of the line pro- 
files, a bright equatorial band results. The algorithm does 



not allow negative image values, which means it has greater 
freedom in setting the intensity (and hence the location) of 
the bright band. The bright band therefore borders the dark 
band, which lies on the equator. 



3.4 Time variability 

One of the key assumptions of our model is that features 
on the surface of the secondary do not vary during the 
course of an observation, e.g. a spot does not suddenly ap- 
pear/disappear at a particular binary phase. Here we ex- 
plore how a variation in the line flux from one region of 
the secondary may effect the reconstruction by simulating 
the presence of a flare. The flare is situated on the equa- 
tor of the leading hemisphere (phase 0.75) and displaced 
slightly towards the back of the star. It is only present be- 
tween phases 0.8 and 0.92 inclusive, has approximately six 
times the intensity of the surrounding non-irradiated region 
and covers 0.7% of the total surface area. 

In the reconstruction (Image D, figure 0) a bright region 
can be seen at phase 0.75. This region is at roughly the 
same longitude but at a much lower latitude than the actual 
position of the flare. This shift in latitude is easily explained; 
due to the orbital inclination, regions near the south pole are 
visible for shorter periods than the equatorial regions. By 
placing a bright region there, the fitting routine is trying to 
account for the limited phase coverage over which the flare is 
observed. However, in doing so the model must smear this 
feature in order to match the observed radial velocities of 
the flare. In addition, bright 'ringing' can be seen, which 
is due to the limited phase coverage over which the flare is 
observed; an ex plana tion of this effect will be given in more 
detail in section 3. IS. 



3.5 Limb darkening 

In section |i| we mentioned that the line profile is scaled to 
take into account limb darkening. The fitting procedure can- 
not, without prior information, account for limb darkening. 
This is because limb darkening causes each element in the 
map to have a phase-dependent intensity variation. We must 
therefore supply the fitting procedure with an appropriate 
limb darkening relation. In this section we show the effects 
of assigning incorrect limb darkening strengths and laws. 

First, data was generated assuming a linear limb dark- 
ening coefficient of 0.5. This was then fitted assuming no 
limb darkening (Image E, figure |l]). The spot features are 
readily visible, although a dark equatorial band is formed. 
This band can be understood because the model calculates 
much higher intensities at the edges of the profiles than are 
present in the fake dataset. As described in section 3.3, this 



results in dark equatorial bands in the reconstruction. In ad- 
dition, the overall intensity level of the map is less than that 
of the test image due to the reduction in light caused by the 
limb darkening. Conversely, if one assumes too high a limb 
darkening coefficient then the opposite occurs. The model 
calculates much lower intensities at the edges of the pro- 
files and in order to compensate the fitting routine produces 
bright bands around the equator. 

For our second test, we assumed an incorrect limb dark- 
ening law. Data generated using a linear limb darkening 
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coefficient of 0.5 was then fitted assuming quadratic limb 
darkening co-efficients of a = -0.025, b = 0.635 and 

I( M )=/i(l-a(l- M )-&(l- M ) 2 ), 

where fi = cos 7, 7 is the angle between the line of sight 
and the emergent flux and I 1 is the monochromatic specific 
intensity at fj,=l. The reconstruction showed a bright band 
around the equatorial regions where the fitting routine has 
had to compensate for the increased limb darkening pre- 
dicted by the quadratic limb darkening law in this case. 



3.6 Velocity Smearing 

Velocity smearing is the degradation in spectral resolution 
due to the combined effects of a finite exposure time and the 
orbital motion of the secondary star. Velocity smearing can 
be ignored if it is insignificant compared to the instrumen- 
tal resolution. To obtain adequate signal-to-noise, however, 
exposure times, i cxp , are usually selected which match the 
resulting velocity smearing to the instrumental resolution, 
v rcs , via the formula 



tcxp ~ Pv ICS /2T:Kft, 



(1) 



where P is the orbital period of the binary and Kr is 
the radial-velocity semi-amplitude of the secondary star. In 
cases where longer exposure times are unavoidable, signif- 
icant velocity smearing will occur. It is possible to correct 
for the effects of velocity smearing in Roche tomography by 
calculating the line profiles at intermediate phases during 
an exposure and averaging the results. This slows the iter- 
ations, however, so in this section we explore what happens 
if velocity smearing is neglected during the reconstruction 
process. 

Image F (figure shows what happens when an ex- 
posure time twice as long as that required by equation [l] 
is assumed when generating fake data and then ignored in 
the reconstruction process. Once again the largest effects are 
seen in the equatorial regions. Figure ^ shows why this oc- 
curs. The radial velocities of the edges of the smeared profile 
extend beyond the velocities expected when the exposure 
time is ignored. For the same reasons as discussed in sec- 
tion 3.3, this results in bright and dark artefacts around the 
equator. Figure |§] also shows how bumps in the line profile 
due to star-spots (e.g. the dip around kms -1 ) are smeared, 
resulting in a blurring of the star-spots in the reconstruction. 



3.7 Intrinsic line profile width 

We have so far assumed that the intrinsic line-profile width 
is negligible in comparison with the rotational broadening. 
This may not always be the case and the width or shape of 
the intrinsic line profile may be important. Here we exam- 
ine the effects of assuming an incorrect intrinsic line-profile 
width. 

Data was produced assuming 20 kms -1 instrumental 
resolution as before but with an intrinsic profile width of 
45 kms -1 (FWHM). The reconstruction was carried out as- 
suming no intrinsic line-profile width and the result was very 
similar to Image F (figure [l]). Once more the edges of the 
lines cannot be fitted and banding is hence seen around the 
equator. The reconstruction also shows that the spots are 




Figure 3. Solid line: line profile at phase 0, computed assum- 
ing zero exposure time. Crossed line: line profile around phase 
0, computed assuming an exposure time twice that required by 
equation jj]. 



broadened. This is due to the algorithm attempting to match 
the width of the bumps in the line profile which, in the ab- 
sence of intrinsic broadening, it can only do by increasing 
the size of the spots in the reconstruction. 



3.8 Eclipse by an accretion disc 

The highest inclination CVs may show signs of an eclipse of 
the secondary star by an accretion disc around phase 0.5. 
This can be accounted for during the reconstruction process 
either by removing the affected spectra (as commonly done 
in Doppler tomography during primary eclipse), or by in- 
cluding the eclipse in the algorithm. If neither of these steps 
are taken, artefacts will be introduced in the reconstructions, 
as described in this section. 

A model including a cylindrical accretion disc of ra- 
dius 0.5Rlx and height O.lRdisc was created and a synthetic 
trailed spectrum including the eclipse by the disc was gen- 
erated assuming i = 80°. Image G (figure [j]) shows the re- 
construction if the eclipse by the accretion disc is ignored. 
As one might expect, the artefacts generated by the eclipse 
are most apparent around the Li point. In particular, the 
southern regions (where light is blocked by the eclipse) ap- 
pear dark, and there are alternating regions of bright and 
dark bands radiating from the Li point. These are due to 
the attempts of the fitting algorithm to match the eclipse 
phases, resulting in a reduction of intensity on the inner 
hemisphere which is compensated for at other phases by 
the bright bands. Note that the spot features are well re- 
constructed, but exhibit the mirroring effect described in 
section 3.11 due to the high inclination. 



3.9 Masses 

A knowledge of the binary masses are essential in order to 
define the geometry of the secondary star model. Very few 
accurate mass determinations are available for CVs, how- 
ever, and though we can use Roche tomography to determine 
them (see section ^|) we still need to know how an error in 
the masses may effect the reconstructions. 

Image H (figure ^) shows the effect of underestimating 
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the secondary mass by O.IM©. The effect of reducing the 
secondary's mass whilst keeping the period and inclination 
the same is to shift the model to higher radial velocities. 
Hence the spot features on the reconstruction appear shifted 
towards the Li point and lower radial velocities in order to 
compensate for this. In addition, as the radial velocities of 
the spot features in the original dataset no longer correspond 
to a well defined location on the secondary, the spots appear 
smeared over a larger area. At phases 0.25 and 0.75 the outer 
hemisphere of the model is shifted to higher radial velocities 
than present in the dataset and hence a dark band is formed 
near the rear of the star. At other phases the model cannot 
match the width of the line profiles due to the reduced radius 
of the star and this results in a bright equatorial band. 



3.10 Inclination 

Up until now we have always assumed the correct orbital 
inclination (60°) in each of our tests. In reality, however, we 
may not know the inclination accurately. Image I (figure [l]) 
shows the effect of assuming an incorrect inclination of 55°. 

Once again, the most obvious artefacts appear at the 
equatorial regions. The lower orbital inclination shifts the 
model to lower radial velocities. As a consequence, the great- 
est mismatch between the observed and computed trailed 
spectra occur around phases 0.25 and 0.75. At the lowest 
radial velocities expected by the model there is no data and 
a similar scenario occurs to that discussed in section 3.3 



dark bands are fitted around the equator on the inner hemi- 
sphere of the star to compensate for the lack of data at these 
velocities. Bright regions then appear around the dark band 
in order to fit the intensities seen at other phases. At higher 
radial velocities around phases 0.25 and 0.75, the observed 
intensity is far greater than expected by the model, and the 
opposite occurs on the outer hemisphere of the star. 

The spot features remain visible, but are systematically 
shifted towards the outer hemisphere of the star. This shifts 
them to regions of the star with higher radial velocities, 
which is necessary to counteract the shift of the model to 
lower radial velocities caused by assuming an incorrect in- 
clination. The opposite occurs if too high a value for the or- 
bital inclination is assumed; the dark bands become bright 
and the spot features are systematically shifted towards the 
inner hemisphere and lower radial velocities. 



also shows that the features become less apparent due to the 
fact that they are smeared over twice the area they originally 
covered. 



3.12 Phase sampling 

Up until now our simulations have been carried out over 50 
evenly sampled phase bins covering a whole binary orbit. In 
practice, however, it may not always be possible to obtain 
this number of phases due to signal-to-noise requirements, or 
there may be gaps in the orbital coverage because of weather 
conditions. Here we model such scenarios. 

First we explore how under-sampling in phase can af- 
fect the reconstruction. Image K (figure [i]) was reconstructed 
from 6 evenly spaced phases (5 independent) using zero ex- 
posure lengthsf] and it is dominated by ring- like streaks. 
The effect is analogous to the stre aks observed in Doppler 
tomography (Marsh & Home 1988) and can be understood 
by considering that lines of constant radial velocity on the 
secondary star at a particular phase can be integrated along 
to construct a line profile. These lines of constant radial 
velocity are ring-like in shape and if there are only a few 
phases, or if the profile is particularl y br ight at a certain 
phase (as in the case of a flare; section 3.4), the streaks will 
not destructively interfere, leaving ring-like artefacts on the 
Roche tomogram. To demonstrate this effect we produced 
maps where those elements with the same radial velocity as 
the spot features were represented as a grey-scale whilst all 
other elements were left blank. This was done for each of the 
6 phases in the dataset and the resulting images were then 
superimposed on top of one another so that regions where 
streaks overlap are shown darker (Image L, figure [j]). It can 
be seen that the features in this image closely resemble the 
artefacts in Image K (figure |l|) . 

Image M (figure gj) shows the effect of incomplete phase 
coverage, with phases between 0.4 and 0.6 missing (which 
may occur if certain phases are discarde d d ue to, for ex- 
ample, an eclipse by the disc; see section 3.8). There is no 
problem recovering the features around the Li point where 
the majority of the information is missing. Although the 
other spot features are also readily recovered they appear 
elongated or streaked due to the same process described in 
the above paragraph. 



3.11 The mirroring effect 

The following sections describe artefacts that are not strictly 
due to systematic errors but are a result of factors that are 
either within our control (e.g. instrumental resolution, phase 
sampling) or are inherent to the reconstruction process itself. 

We begin with the problem that, although the radial 
velocities contain latitudinal information, they cannot con- 
strain whether a feature is located in the northern or south- 
ern hemisphere. This information can only be supplied by 
self-obscuration, but this becomes ambiguous at higher in- 
clinations until, at an inclination of 90°, a feature in the 
northern hemisphere is self-obscured in an identical manner 
to a feature located at the same latitude in the southern 
hemisphere. This results in a 'mirroring' of features in both 
hemispheres as demonstrated in Image J (figure [l]), which 



3.13 Instrumental resolution 

The spectral resolution of the data governs the size of the 
features that can be imaged, as shown by the blurred spots 
in Image N (figure jj]) where we have reconstructed the test 
image using data of lower resolution (60 kms -1 ) than be- 
fore. The number of surface elements in the model should 
be chosen to match the spectral resolution. If the model has 
too few elements, the separation between adjacent elements 
in velocity space (particularly around phases 0.25 and 0.75) 
is greater than the resolution of the data. This means that 

* Although such a scenario would probably only occur if long 
exposure times were being used, we have assumed zero exposure 
length in order t o sep arate the effects of velocity smearing (de- 
scribed in section |3.6|) and those due to phase under-sampling. 
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there will be parts of the line profile to which no surface ele- 
ments in the model contribute (at a particular phase) and no 
acceptable fit will be found. If, on the other-hand, too many 
elements in the model are used, there will be no problem in 
fitting the data but it will slow the reconstruction process. 



3.14 Default maps 

In the reconstruction process the final map is the one of max- 
imum entropy (relative to an assumed default map) which 
is consistent with the data. If the data are good then the fi- 
nal map is constrained by the data and will not be strongly 
influenced by the default. In this case the choice of default 
makes little difference to the final map. If the data are noisy, 
however, the data constraints will be weak and the map will 
be strongly influenced by the default. The default can there- 
fore be regarded as containing prior information about the 
map, and regions of the map that are not constrained by the 
data will converge to the default value. 

To investigate the effect of different default maps on the 
reconstruction, it is necessary to degrade the quality of the 
test dataset. We therefore generated data over 15 phase bins 
with a resolution of only 60 kms~ and a S/N ratio of only 
50. The dataset was then fitted using a uniform default map 
and a Gaussian-blurr default map. The former was created 
by setting every element in the default to the average value 
of the map, and hence constrains large-scale surface struc- 
ture. The latter was created by smoothing the map using a 
Gaussian-blurring function and hence constrains small-scale 
surface structure. 

Image O (figure |l|) shows the reconstruction using the 
uniform default map and, as expected, it can be seen that 
the spots are more blurred than in the reconstruction using 
the Gaussian-blurr default map (Image P, figure |l]). Fea- 
tures below -60° latitude are never visible and hence in 
the fitting routine they are assigned a value determined by 
the default map. In the case of the uniform default map, 
this is the average element value of the map. In the case of 
the Gaussian-blurr default map, however, the value assigned 
depends upon the intensity values of the neighbouring ele- 
ments and, for elements that are not visible, this is close to 
zero intensity. This explains why the reconstruction using 
the Gaussian-blurr default has a dark southern limb and, 
to compensate for this, a bright region above the Li point. 
Hence, although the use of a Gaussian-blurr default is more 
likely to resolve small features, it is important to be aware 
of the effects that visibility may have on the final outcome. 



4 DETERMINING PHYSICAL PARAMETERS: 
THE ENTROPY LANDSCAPE 

If the centre-of-light and centre-of-mass of the secondary 
are not coincident, the star's radial- velocity curve will be 
distorted in some way from the pure sine wave which repre- 
sents the motion of its centre-of-mass. The observed radial- 
velocity curves of CV secondaries suffer from this distortion, 
due to both geometrical effects caused by the Roche-lobe 
shape and non-uniformities in the surface distribution of the 
line strength due to, for example, irradiation. If these fac- 
tors are ignored, they will lead to systematic errors in the 
determination of the binary star masses. 



By modelling the shape and mapping the surface inten- 
sity distribution of the secondary star, Roche tomography 
can provide constraints on the CV masses. This is because, 
as discussed in section 3.9, incorrect values of Mi and Mi 



will introduce artefacts in the reconstructions which will 
lower the entropy of the final solution. The correct values 
of the component masses are therefore those that produce 
the map of highest entropy. 

This is demonstrated in figure ^j, which displays the 
'entropy landscape' of the test data shown in figure H. Each 
square in the entropy landscape corresponds to the maxi- 
mum entropy obtained in a reconstruction assuming a par- 
ticular pair of component mass values, a uniform default 
map and iterating to the same x 2 value on each occasion. 
The reconstruction with the highest entropy was obtained 
for component masses of Mi = 1.0 Mq and M2 = 0.5 Mq, 
which are identical to the masses used to construct the test 
data. This shows the entropy-landscape method to be an 
extremely effective method of determining accurate compo- 
nent masses which accounts for both geometrical distortion 
and complicated intensity distributions on the secondary 
star. Note that the correct inclination was used for every re- 
construction in figure but in practice the inclination may 
not be accurately known. A similar method to the entropy 
landscape can then be used to determine the inclination or, 
in cases where the secondary eclipses the white dwarf, the 
inclination can be varied for each combination of compo- 
nent masses in order to remain consistent with the observed 
eclipse width. 

Figure [l| also shows a ridge of high entropy running 
along a line of (almost) constant q. Component masses lying 
along this diagonal are therefore more difficult to distinguish 
between than component masses lying perpendicular to it. 
The gradient of the diagonal can be derived by considering 
that the best fit is obtained when the centre of the model 
line profile matches the centre of the observed profile. The 
centre of the line profile is governed by the radial velocity 
of the secondary, which relates directly to the distance of 
the centre-of-mass of the secondary from the centre-of-mass 
of the binary, 02. Thus only those combinations of Mi and 
M2 that maintain constant 02 can provide satisfactory fits 
to the data. Using Kepler's law and assuming a constant 
period and at we can show that the ridges in the entropy 
landscape correspond to component masses which obey the 
relation 



Mi 



= constant. 



(2) 



(Mi + M 2 ) 2/3 

This si mple argument is complica ted in the case of irra- 
diation. Rutten & Dhillon (1994) found that the ridge of 



high entropy was shifted towards slightly higher values of 
M2. Their explanation for this was that the intensity distri- 
bution was strongly peaked towards the Li point, causing 
a large difference with the uniform default map. This dif- 
ference was reduced by increasing M2, allowing the bright 
region to shift away from the Li point and to spread over 
more elements, reducing the deviation between the map and 
the default whilst still fitting the trailed spectrum to within 
its uncertainties. Despite this, the mass determination was 
still correct to within ~2 per cent. 

Another feature of note in the entropy landscape of fig- 
ure ^ is that, as the primary mass increases, there is a larger 
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Figure 4. Entropy landscape of the test data presented in fig- 
ure ^. The contours are lines of constant entropy and have been 
plotted to highlight the ridge of high entropy. 



range of secondary masses that can be fitted satisfactorily. 
This can be explained by generalising equation |^ to 

g+lccM^ 2 , 

which shows that as the primary mass increases, the mass 
ratio, q = M2/M1, must also increase in order to hold a 2 
constant. It can be shown that the secondary star radius is 
related to q and 02 via the equation 



R2/C12 <xq 1/3 (l + 



,2/3 



which shows that the radius of the secondary star must in- 
crease if q increases. This in turn means that the model 
line-width is greater than the observed line-width and this 
overlap allows a greater range of secondary star geometries 
for which a satisfactory fit can be obtained. In contrast, at 
low primary masses the allowed value of q is decreased, re- 
sulting in a low secondary star radius and hence a small 
model line-width. At the lowest primary masses in figure 
the model line- widths become too narrow to fit the observed 
profile to the desired x 2 > no matter what intensity pattern 
is mapped onto the secondary star's surface, explaining the 
abrupt end of the high entropy ridge on the lower left-hand 
portion of figure ^. 

It is important to note that the maximum entropy al- 
gorithm we use requires that all input data is positive. As a 
result, it is necessary to invert absorption line profiles and, 
in the case of low signal-to-noise data, add a positive con- 
stant in order to ensure that there are no negative values 
in the continuum-subtracted spectra. If this positive con- 
stant is not added, and any negative points are either set 
to zero or ignored during the reconstruction, the fit will be 
positively biased in the line wings, resulting in a broader 
computed profile and hence larger stellar masses in the en- 
tropy landscape. It is a simple matter to account for this 
positive constant during the reconstruction process by em- 
ploying a virtual image element which contributes a single 
value to all data points, effectively cancelling out the con- 
stant. This virtual image element can also be used to correct 
for small errors in the continuum subtraction, which might 
otherwise systematically affect the masses derived from the 
entropy landscape. 



5 STATISTICAL ERRORS 

Although it is possible to propagate the statistical errors on 
the data points through the iterative process in order to cal- 
culate the statistical errors on each element in the map, the 
resulting error-bars are unreliable due to the fact that noise 
in Roche tomograms is correlated. This correlation arises 
from the projection of bumps in the profile along arcs of 
constant radial velocity on the secondary star and from the 
effects of blurring by the default map. Doppler imaging of 
single stars has largely relied on consistency tests, such as 
comparing ima ge reconstructions using either subsets of ob- 
serva tions (e.g. Barnes et al. 1998 ) or different spectral lines 
(e.g. Strassmeier & Bartus 2000), in order to address this 



issue. So far, the only attempt at a formal statistical error 
analysis in Doppler i maging has been u ndertaken using the 
Occamian approach ( Bcrdyugina 1998| ). 

In the case of Roche tomography, one will undoubtedly 
be working with lower signal-to-noise datasets and fewer 
spectral lines than conventional Doppler imaging. For in- 
stance, the well-st udied K-star AB Dor h as a mean visual 
magnitude of ~ 7 (Cameron & Foing 1997), whereas the sec- 
ondary stars in CVs never attain magnitudes of below ~10 
and are usually several magnitudes fainter than this. Roche 
tomograms of CV secondaries are therefore not as strongly 
constrained as Doppler images of single stars by the input 
data. This makes the determination of statistical errors of 
prime importance in Roche tomography, because only then 
can one test whether a surface feature is real or just a spu- 
rious artefact due to noise. 

We have addressed the problem of statistical error de- 
termination in Roche tomography using a Monte-Carlo style 
approach. Monte-Carlo techniques rely on the construction 
of a large (typically hundreds, in the case of Roche tomogra- 
phy) sample of synthesized datasets which have been effec- 
tively drawn from the same parent population as the original 
dataset, i.e. as if the observations have been repeated many 
hundreds of times. This large sample of synthesized datasets 
is then used to create a large sample of Roche tomograms, 
resulting in a probability distribution for each element in the 
map. The main difficulty with this technique, aside from the 
demands on computer time, lies in the construction o f the 



sample of synt hesized datasets. One approach (used by Rut 



ten et al. 1994 in spectral eclipse mapping) is to 'jiggle' each 
data point about its observed value, by an amount given by 
its error bar multiplied by a number output by a Gaussian 
random-number generator with zero mean and unit vari- 
ance. This process adds noise to the data, however, which 
means that the synthesized datasets are not being drawn 
from the same parent population as the observed dataset 
- noise is being added to a dataset which has already had 
noise added to it during the measurement process. In prac- 
tice, this means that fitting the sample datasets to the same 
level of x 2 as the original dataset is either impossible (i.e. 
the iteration does not converge) or results in maps domi- 
nated by noise, which grossly overestimates the true error 
on each element. We have verified that this is the case using 
simulations similar to those described in section [0 . 

Instead we hav e opted for the bootstrappi ng technique 



of |Efron (1979)| see |Efron fc Tibshirani (1993)1 for a general 



introduction to the method. We implemented the bootstrap 
as follows. From our observed trailed spectrum containing n 



© 0000 RAS, MNRAS 000, 000-000 



10 C.A. Watson and V. S. Dhillon 



data points we form our synthesized trailed spectrum by se- 
lecting, at random and with replacement, n data values and 
placing these at their original positions in the new synthe- 
sized trailed spectrum. For points that are not selected we 
set the associated error bar to infinity and thereby effectively 
omit these points from the fit. For points that have been se- 
lected once or more than once (typically 37% of the points) 
we divide their error bars by the square root of the number 
of times they were picked. The advantage of bootstrap re- 
sampling over jiggling is that the data is not made noisier 
by the process as only the errors bars on the data points are 
manipulated (i.e. the data values remain unchanged) . In so 
doing, the amount of noise present in the data is conserved 
and it is therefore possible to fit the synthesized trailed spec- 
tra to the same level of x 2 as the observed data, giving a 
much more reliable estimate of the statistical errors in the 
maps. 

5.1 Testing the bootstrap 

In order to assess the viability of the bootstrap technique 
we ran a simple test. A model with parameters Mi = 1.0 
M , M 2 = 0.5 M Q , i = 60° and P = 2.78 hours was cre- 
ated of uniform intensity distribution, with the exception of 
a single spot covering 0.8% of the total surface area on the 
trailing hemisphere (see figure ^). From this model a syn- 
thetic trailed spectrum was generated with 30 evenly spaced 
phase bins and 60 kms" 1 instrumental resolution. Each data 
point in the spectrum was then 'jiggled' by 10% of the orig- 
inal data value in order to simulate noise. A close fit to this 
dataset was then carried out using the correct parameters. 

As can be seen, the reconstruction (figure ^j) shows 
many spurious artefacts due to noise. Fitting to a higher 
X 2 , whilst removing some of the lower- level noise artefacts, 
does not remove the spurious spot features marked A and 
B in figure []. The triangular symbols in figure ^ show the 
intensities in each vertical slice through A, B and the real 
spot. It can be seen that, from this information alone, it is 
impossible to deduce which spot features are real. 

The curves in figure (7] show the confidence intervals 
found after reconstructing 200 bootstrap samples of the im- 
age. As the distribution of element intensities found is often 
non-normal and not centred on the intensity found in the 
original fit (see figure ^ , it is incorrect to calculate a sum- 
mary error statistic like the standard deviation. Instead we 
take the mode of the distribution to represent the 'most 
probable' value, and define the 95% confidence interval (for 
example) as the region which encloses 95% of the bootstrap 
values above and below the mode. It can be seen that the 
confidence intervals in figure |^ widen significantly around 
the noise features A and B, and encompass the true inten- 
sity level (indicated by the horizontal line). The confidence 
intervals around the real spot feature, however, do not ex- 
hibit such a widening and clearly do not encompass the true 
intensity level of the non-spotted regions. This test there- 
fore confirms that the bootstrapping method works; it has 
correctly identified the real spot feature and correctly iden- 
tified features A and B as artefacts due to noise. In a similar 
test, statistical errors determined by the method of 'jiggling' 
failed to distinguish between the real spot and features A 
and B. 

There are two additional features of note in figure Q 



= 0.25 




= 0.75 




North Pole 



saturated 




Figure 5. The test im age used to construct the noisy dataset 
described in section 



= 0.25 



= 0.75 



= 0.5 



North Pole 



saturated 





Figure 6. The Roche tomography reconstruction of the test im- 
age in figure ^| using noisy data. Although the spot can be seen at 
phase 0.25, many spurious artefacts due to noise are also present, 
including features A and B. 



First, it can be seen that in regions where the data do not 
constrain the intensity distribution (such as the southern 
hemisphere, which is never visible), the reconstruction al- 
ways converges to the default map value, which is why the 
confidence intervals converge around 'SP' in figure (?]. Sec- 
ond, the slices through features A and B demonstrate how 
noise is correlated in the Roche tomography reconstruction 
process, because elements surrounding the peak intensity are 
all shifted in the same direction as the peak itself. 



6 CONCLUSIONS 

We have shown that any feature on a Roche tomogram 
must be subjected to two tests before its reality can be con- 
firmed. The first test is to determine whether the feature 
is statistically significant and is performed using a Monte- 
Carlo technique. The conventional method, where hypothet- 
ical datasets are constructed by 'jiggling' the data in accor- 
dance with their error bars, was found to add noise and 
hence grossly over-estimate the errors in the reconstruction. 
The solution to this problem was found by using a bootstrap 
re-sampling algorithm and, through the use of simulations, 
we have shown this technique to correctly distinguish be- 
tween real features and artefacts due to noise. The second 
test is to compare the feature with the appearance of known 
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Figure 7. Triangles: intensity values along vertical slices through 
A (top), B (middle) and the real spot (bottom) in figure^. Curves: 
confidence intervals along vertical slices through A (top), B (mid- 
dle) and the real spot (bottom) in figure ^. LH, NP, TH and SP 
at the top of the figure represent the positions of elements at the 
leading hemisphere, north pole, trailing hemisphere and south 
pole, respectively. 



intensity (arbitrary units) 



Figure 8. Histogram of intensities obtained from 200 boot- 
strapped reconstructions for one element. The vertical line rep- 
resents the intensity value obtained in the original fit. The true 
value is 61. 



artefacts of the technique. In general, we find that system- 
atic errors result in only three major artefacts in the Roche 
tomograms: ring-like streaks, equatorial banding and blur- 
ring. The effects of these artefacts can be minimised, how- 
ever, by fine tuning the input parameters. A fine example 
of this is given by the entropy landscape technique, which 
also provides accurate component masses. If a surface fea- 
ture survives both of the above tests unscathed, it can be 
assumed to be real. 
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